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Abstract. A time evolution operator in the interaction picture is given by exponentiat¬ 
ing an interaction Hamiltonian H. Important examples of Hamiltonians, often encoun¬ 
tered in quantum optics, condensed matter and high energy physics, are of a general 
form H = r{A' — A), where A is a multimode boson operator and r is the coupling 
constant. If no simple factorization formula for the evolution operator exists, the calcu¬ 
lation of the evolution operator is a notoriously difficult problem. In this case the only 
available option may be to Taylor expand the operator in r and act on a state of interest 
1 ^. But this brute-force method quickly hits the complexity barrier since the number 
of evaluated expressions increases exponentially. We relate a combinatorial structure 
called Dyck paths to the action of a boson word (monomial) and a large class of mono¬ 
mial sums on a quantum state This allows us to cross the exponential gap and 
make the problem of a boson unitary operator evaluation computationally tractable by 
achieving polynomial-time complexity for an extensive family of physically interesting 
multimode Hamiltonians. We further test our method on a cubic boson Hamiltonian 
whose Taylor series is known to diverge for all nonzero values of the coupling constant 
and an analytic continuation via a Fade approximant must be performed. 


1. Introduction 

The calculation of an evolution unitary operator V = exp [—tH], where H is an in¬ 
teraction Hamiltonian in the interaction picture, is a generic problem encountered in 
quantum field theory, quantum optics or condensed matter physics, to name a few. Of 
considerable interest are single and mainly multimode boson Hamiltonians describing 
the interaction between two and more boson field modes. The Hamiltonian H is a Her- 
mitian operator and can be written as H = vA' + vA, where v e C and the bar denotes 
complex conjugation. The operator A is a function of one or more boson creation and 
annihilation operators Uj and a.', respectively, satisfying the Weyl-Heisenberg algebra 
[uj, aj ] = 5ij id. The index i denotes the i-th boson mode and id is an identity operator 
(not necessarily represented by a finite-dimensional identity matrix) commuting with U; 
and a.. Because the Hamiltonian is often a sum of two or more non-commuting terms, 
the problem of finding an explicit form of V is not a straightforward one. The standard 
way of solving it is to use perturbation theory by expanding V as the Dyson series [|li [^. 
It can also can be converted to the problem of factorization of the operator exponen¬ 
tial function knovm as the BCH formula (after Baker-Campbell-Hausdorff) [j^. Strictly 
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speaking, it is the Zassenhaus formula [Q, [H] that provides an explicit way of factorizing 
operator or matrix exponentials but it is not the only, see [ji -11 ]. This extensive topic 
is also reviewed in [12] and more recently in [f^fl^ . 

By setting v = s + ir we will restrict our attention to the evolution operator 


V = exp [r(A' — A)] . 


( 1 ) 


From all the existing approaches to evaluate Eq. dl]), the Taylor series expansion around 
r = 0 happens to be (most likely) the least sophisticated attempt to compute the ex¬ 
ponential. But if no simple factorization theorem exists (for example when A and A' 
commute) and other methods fail or are equally inefficient, it may be the only available 
option. As a matter of fact, short-time evolution is often studied by calculating the first 
few terms of a Taylor expansion. However, for longer times the exponential growth in 
k of the number of summands (A' ± A) quickly makes such a calculation intractable 
(note that we do not claim that all methods to evaluate V are exponential in time). 

The situation could have been saved by the following procedure. The expression 
(A' ± A)^ can be brought to normal form, where all creation operators aj, forming A 
and A^, are to the left of all annihilation operators U;. We will denote the corresponding 
operation by Af. The main advantage of normally ordered operators lies in an easy 
calculation of their action on many states of interest, in particular, the eigenvector of a 
(a coherent state |a) with the Minkowski vacuum as its special case). However, in 
general it is not easy to find the normal form. As an example, take a boson word (also 
called product, string or monomial) w = a^a^^a. The action of J\f results in 

M[w} = + 6a'a. (2) 


A straightforward, but not recommendable, way to obtain ([2]) is to commute through 
the forest of operators. Eventually, one can use shortcuts in the form of derived com¬ 
mutation rules to make the process less painful |]2[] . A better way is to take advantage 
of a clever differential representation of the aforementioned algebra [Q] 

5 

a a + —, a' ^ a, (3) 

da 


where the bar denotes complex conjugation. Then the above calculation can be done in 
a blink of an eye (the celebrated Eock-Bargmann representation is another option [1141 ) 


but this is the exception rather than the rule. In the case of A/^[(A^ ±A)^], the com¬ 
mutation relations, the differential operator representation ([3]) or any other technique, 
such as Wick’s theorem i] , would lead to the right answer for but it may be difficult to 
obtain a general form for an arbitrary A and k or the computational complexity of the 
J\f operator can be insurmountable. 

The problem of boson normal ordering has a long history and a number of techniques 


were developed [15-2^]. Certain special cases of A/” [(A'±A)^] were previously studied 
resulting in elegant combinatorial formulas. But the general treatment for multimode 
boson operators is missing and the calculation of evolution operators generated by such 
Hamiltonians is the main goal of this paper. 


The idea behind our technique is to study (A' ±Aj , and consequently the expansion 
of Eq. ([1]), in the weak sense, that is, by acting on a state of interest. After all, this is the 
situation most often encountered when solving physical problems - the state is usually 
a vacuum state jvac) (or a ground state defined in general by Ajvac) = 0) and one is 
interested either in the vacuum expectation value of an operator or a unitary evolution 
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of the state it acts on. So our method lacks the generality of the factorization theorems 
for Eq. ([!]) as a standalone operator [Isils, 13] or of the formulas developed to normal 


order certain special cases of (A' ±A) 
enables us to analytically calculate 


However, the virtue of our technique is that it 


(A'^ ±A) 


( 0 ) 


(4) 


for an arbitrary multimode boson monomial A, where — |vac) defined above as 
A|vac) = 0 or any state from the semi-infinite or finite tower of states generated by 
the repeated action of A' on |vac). In this work we will focus on the vacuum state 
as it the most often encountered scenario. Equally importantly, the calculation will 
be shovra to be tractable since it can be executed with polynomial-time complexity 
Oidk^) if A' or A is a monomial of the length d. The method does not depend on 
a particular algebraic structure of A apart from the one inherited from a, a^ satisfying 
the Weyl-Heisenberg algebra. Our approach is combinatorial. We will show that there 
exists a common structure called a generalized Dyck path, which is a subject of study in 
analytic combinatorics, computer science and stochastic processes, among others |]25l - 
[ 2 ^ . This insight will lead us directly to a recursive summing formula allowing to 
analytically calculate Eq. ([4|) without the need for normal ordering. The most important 
consequence is that it can be used to efficiently find (in polynomial time) an analytical 
expression for the expanded evolution operator ([1]) 





for any choice of 0 < fC < 00 and 0 < r < 00 where again = |vac) or a tower of 
states generated by A' acting on |vac). 

Multiboson Hamiltonians studied in this paper were often identified as quasi-exactly 


solvable [^IsoD with a close relation to the boson realizations of polynomial alge¬ 
bras [13 iL 13 211 . Among the applicable physical scenarios, Bose-Einstein condensates and 


nonlinear optical systems are the most prominent. Quasi-exactly solvable models have 
been extensively studied [33-35[] and the difference compared to exactly solvable mod¬ 
els is that the spectrum problem of such a Hamiltonian cannot be formulated as the 
usual eigenvalue problem in matrix algebra. However, there exists an invariant sub¬ 
space of functions that are mapped to themselves and so it can be “decoupled” from the 
rest of an otherwise infinite-dimensional matrix Hamiltonian. This finite subspace can 
subsequently be analytically diagonalized and part of the spectrum is recovered. 

Generalized Dyck paths are introduced in Sec. 12.1! and some of their combinatorial 
properties are summarized. In Sec. 12.21 we link Dyck paths to the action of boson 
words, show that all monomials behave as ladder operators and uncover a large class 
of sums of multimode boson monomials that behave as ladder operators as well. This 
enables us to derive our main result in Sec. 12.31 the evaluation of Eq. ([4|) for monomials 
and their sums. In Sec. |2.4l we discuss the computational complexity of our approach. 
We illustrate the evaluation on several examples of boson operators A and compute 
Eq. ([4|) for low values of k (to be able to verify the results by hand) and for some 
very high values of k (to show the speed-up offered by our method). We also show 
that the fast evaluation of Eq. Q leads to an efficient calculation of a boson evolution 
operator for previously hard-to-reach values of the coupling constant r if calculated by 
a series expansion. Our final example is the evolution operator for the three-photon 
degenerate process whose Taylor expansion is known to diverge for any value of the 
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Figure 1. We plot (a) Dyck path D(8,0,0) (b) more general Dyck path 
D(10,0,3) and (c) the generalized Dyck path D(7,3,2). Dyck paths are often 
described as mountain ranges. 


coupling constant. Yet our approach is still useful for the calculation of an analytically 
continued unitary operator using the Fade approximation proposed in . 


2. Generalized Dyck paths and their combinatorics 


2.1. Counting Dyck paths. Let’s introduce an integer lattice path known as a Dyck 
path and some of its properties [251-28]. The lattice for Dyck paths is the set of all 


integer points N x N c A Dyck path D(k) is a path starting at (0,0) and ending at 
(fc, 0), such that the only allowed steps (also called segments or letters in this paper) 
are U = (1,1) and D — (1,-1), that is, going to the nearest integer point northeast or 
southeast. The constraint on a Dyck path is that it never falls below the x axis, but it 
may touch it at any number of integer points between 0 and k. A more general Dyck 
path T){k, 5) starts at (0,0) but is allowed to end at any point (k, 5) following the same 
constraint. Finally, what we will call a generalized Dyck path D(fc, 5^,52) is a lattice 
path starting at (0,5^) and ending at (k, ^ 2 ), where k, 81,82 > 0, and as before, the 
path is constrained to the positive quadrant including the x axis. The value of 5^ can 
be any non-negative integer but the value of 52 depends on k and 5^. Assuming 5^ 
even, if k is even (odd) then only even (odd) values of 82 are possible. For 5^ odd, if k 
is even (odd) then only odd (even) values of 52 are possible. 

Two examples of a generalized Dyck path are in Fig.[TJ Note that 1)(k, 5) = 2)(k, 0, ^ 
and from now on whenever we say a Dyck path we will mean a generalized Dyck patfo. 
A Dyck path is unambiguously defined by its starting point 5^ and a set of instructions 
called a Dyck word that leads the path through the lattice. It is a string of the segments 
U and D whose length is k and our convention will be to read the string from the right. 
The highest Dyck path D(k, 5^, 82 ) is defined by the (highest) Dyck word, where all the 
U segments are on the right: 


W ^D. 


.DU. 


U. 


k+5i-52 
2 


k-5i+52 
2 


Then, all Dyck paths for the fixed parameters k, 5^ and 82 lie “between” the highest 
Dyck path and the x axis. Equivalently, the lowest Dyck path ®(k, 5^, 82 ) is defined 
such that there is no other Dyck path between this path and the x axis for the chosen 


1 


Sometimes, a Dyck 
and D = (1,-1) itHlsj 


lath is called generalized if two steps are allowed - not necessarily being U = (1,1) 
]. This is not the kind of generalization considered in this paper. 
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Figure 2. All G(6, 0,2) = 9 Dyck paths of length k = 6 starting at 5^ = 0 and 
ending at 52 = 2 are depicted. 


set of parameters. It is generated by the lowest Dyck word 

W = U . UDU . DUD . D. 


fc-5i -5, 


Any Dyck path with the given parameters k, 5^, 52 contains the same number of U and 
D segments. 

The number of Dyck paths T>(k, 0,0) is famously counted by the Catalan number Cfc /2 


|D(k,0,0)| = Cfc/2- 


2 + fc 


(5) 


Catalan numbers are ubiquitous in computer science and discrete mathematics [1391 an d 
this makes Dyck paths a well studied object appearing in many reincarnation^ 250. 
One of them is Bertrand’s ballot problem and by virtue of this insight one finds [l40l1 the 
number of Dyck paths D(k, 5^, 52 ) to be 


im5i,52)| = G(k,5i,52)- 


Uk + 52-5,)J U(/c-52-5i-2);’ 


( 6 ) 


where the binomial coefficient is defined to be zero whenever the “denominator” is 
negative. 

As an illustration we plotted all Dyck paths I>(6,0,2) in Fig.[2]for which the expres¬ 
sion from Eq. @ simplifies to 


G{k, 0 , 82 ) = 


52 + 1 

52+k j 


Sn+k 


(7) 


'A huge collection of examples is available at http://ww-math.mit.edu/~rstctn/ec as an ad¬ 


dendum to 
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Note that if /c = 5^ + ^2 in ® , all the Dyck paths D(k, k — 82 , 82 ) are confined to a 
rectangle whose two corners are touching the positive x and y axes and 

G{k,k- 82 , 82 ) 

is a binomial number. This is nothing else than a binomial walk on the Pascal triangle 
(rectangle) experimentally realized by the Galton board. 


2.2. Dyck paths and bosons. Let A be a general multimode boson operator such that 
there exists xp satisfying A|vac) = 0. So we may call |vac) a vacuum or ground state. 
The object of our study is the expression A' ± A that can be considered as a hypothetical 
Hamiltoniar0. The operator A' acts repeatedly on the ground state and generates a 
finite or eventually semi-infinite tower of states: A'^ |vac) = where |vac) = 
and is, in general, an unnormalized state. The first crucial question is under what 
conditions the adjoint operator A behaves as an annihilation operator: 

A'xp^P~^^ — ?ipXp^P\ ( 8 a) 

Axp^P^ = ^pXp^P~^\ ( 8 b) 

We adopted the convention that the states with no tildes are normalized and so the 
coefficients Xp and are determined by this condition. The second crucial point we 
have to address in order for our method to work is that our approach hinges on the 
ability to analytically sum over the product later derived in Sec. 12.31 and used in 
Eqs. (I39D and (I49D (appearing as X„^.ijl„^.). First, we will show that Xpfip can be obtained 
even without knowing the normalization of xp^P^ (note that we are not interested in Xp 
and jjLp but only in their product). Finding the normalization may not only be rather 
laborious for a generic boson expression A and all p but we even cannot guarantee the 
existence of a closed, or even simple, form. So assume for a moment that A, A^ behave 
as proper ladder operators as in ® . Then we may write 

AA'xp^P~^^ = XpAxp^P^ = XpPLpXp^P~^\ (9) 


But this clearly does not depend on whether xp^P~^^ is normalized or not. It holds that 
. 0 (p-i) _ where A'^ ^ |vac) = xp^P~^^ and is the normalization of 

Hence we get the sought product XpPp from the unnormalized ladder states for 

all p 

AA'i''0(P-i) = XpPpXp^P-^'^ (10) 


greatly simplifying its derivation. 

Let’s now uncover a large class of the boson expressions A, A' that behaves as ladder 
operators. Remarkably, we will show that this is equivalent to the condition that the 
product XpPp is an integer and therefore leads to the announced speed-up derived in 
Sec. [231 


Definition. A normally ordered monomial of the length d is the expression 


= a 


■tfko j 


aiA..a, a/. 


^0 V- 

We define a normal disjoint monomial of the same order as any normally ordered mono¬ 
mial that is created by taking the creation and annihilation exponents (kg, ..., fc^j) 


^The minus version needed to evaluate O is recovered by putting a minus to all products with an odd 
number of A’s. The appearance of the plus sign is irrelevant for the presented method. 
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and and (independently) permuting them such that no mode with a 

nonzero exponent in is present in the newly created one. 

The definition seems restrictive but the opposite is true. It allows for a large num¬ 
ber of multimode boson expressions that we will appreciate in the example after the 
following statement. 

Proposition. Let Abe a sum of disjoint multimode boson monomials of the same order and 
I vac) a vacuum state such that A | vac) = 0. Then the operator AA' is proportional to an 

identity map for all ladder states = A ^ |vac)j that is, AA^ip^P 

Moreover, the constant of proportionality qp-i is an integer. 

Remark. Note that a vacuum (or ground state) is in general not unique. A given oper¬ 
ator A can annihilate several distinct vacui. 


Example 1. The boson expressions A can now be constructed in a great variety. The 
trivial example is every monomial M (arbitrarily ordered). To illustrate the definition 
we further restrict to normal monomials but later we will comment on differently 
ordered alternatives. If we fix the monomial order d then the first nontrivial example is 

= a^ ° for any hg- The disjoint monomials of the same order are all = a! ”, 
where 1 < i < d and so the boson expression A can be any sum 



!=0 


( 12 ) 


Another example is = a^aQa!^^. A sum of disjoint multimode boson monomials of 
the same order can be, for instance, 

-•-3 ---3 ---3 

A = aQaoa| - 1-020203 - 1-040403 . (13) 


To prove the proposition we will need an auxiliary lemma. 


Lemma. Let Mi and M 2 be boson multimode monomials such that M 2 is an arbitrarily 
reordered version of Ml. Then 

Ml -M 2 +/(/!;), (14) 

where /(o;) is a multivariate polynomial with the variables being the i-th mode boson 
numbers n, = a.aifor 0 <i < d. 

Proof Any monomial Mi (normal ordered or not) can always be recast to the form of 
M 2 resulting in Eq. 014D . Then, all summands of the (non-unique) polynomial f(nf) 
must contain an equal number of creation and annihilation operators for all modes 
involved in Mi and M 2 . If it was not the case then the lemma assumption about Mi and 
M 2 would be violated. The polynomial /(o;) can be put (by normal ordering) to a form 
containing only the expressions a!^af and every such building block can be retvritten 
as a sum of powers of O;. ■ 


Proof of proposition. We will be concerned with the sums of normally ordered disjoint 
multimode boson monomials denoted A and A' and their corresponding products AA' 
and A'A. Each monomial in AA' has a reordered counterpart with the same number 
of creation and annihilation operators in A'A. According to the previous lemma, this 
implies 


AA'' =A^A + g(^ni) 


(15) 
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where gCn;) is a multivariate polynomial with the variables being the boson number 
operators Uj. Furthermore, using the distributivity of the commutator and the crucial 
fact that 

I 

and its adjoint A' are disjoint (/ is a disjoint multi-index, see the examples in Eqs. (I12D 
and (I13D ~). we may write 

I 

where [A/,A'|] = f(n[). This is the only surviving term in commutator Eq. (1151) . To use 
this fact, we write 

= AA'^ |vac) (16a) 

= (^A'A +g(nj)j A'^ ^ |vac) (16b) 

= ^A'AA'^ ^ -l-A'^ |vac), (16c) 

where in the second summand of the third row we used the following consequence of 
the boson commutation rule: 

riiaf + (17) 

riia^j (18) 

So we “propagated” the polynomial g(n;) to the right side of A'^ ^ converting it to 
another polynomial g^^~^^(nj), where only the boson number operators appear as vari¬ 
ables. Recall that A' is a sum of disjoint monomials and so all but one /(«/) commutes 
and the one that does not commute switches the side according to the rules Eqs. 017D in 
the same way for all firij). This results in g(n;)Al = A'g*^^^(n;) implying g(n;)A'^ ^ = 
^g(p-i)(nj). We continue by writing 

A'i'AA^^”^ - A'' (aU -F g(n;)) (19a) 

- A''^AA'^^“^ -FA'^''“^g^^'“2^(ni) (19b) 


until we arrive at 


p-i 


AA^^P-^^= A^^A+A^^ \vac) ^ qp_^xp^P-^\ 


( 20 ) 


S = 1 


where we used A|vac) = 0 => A'^A|vac) = 0 and the fact that only the number coeffi¬ 
cients of g^P~^\n j) survive the encounter with |vac) generating qp ■ 


Remark. The A expressions are normally ordered but we might have equally used anti- 
normally ordered monomials [jH, ITsIl or arbitrarily ordered ones and prove the proposi¬ 
tion for them. It would lead to an ever larger class of operators satisfying the proposi¬ 
tion. 


Remark. The class of operators uncovered in the proposition is certainly not the most 
general class of boson expressions with this kind of ladder-like behavior. Yet, among 
this class, the monomials represent a huge variety of hypothetical Hamiltonians and so 
we will mainly focus on them. 
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Example 2. As our first example, consider A' = The ground state is = 

|0)o |n)^ since = 0 for all n and A' (acting repeatedly on generates a finite 

tower of states for 0 < p < n. The corresponding Hamiltonian Hgg — r{^a^ai — 
aouj) is incidentally an important one. It represents one of the generators of the su(2) 
algebra describing the action of a quantum-optical beam-splitter and the fact that the 
ground state is somewhat non-unique turns out to be important too. The I = n + 1 
states generated by A^^ |0)o |n)i carry the Z-th representation of su(2) [14 111 . The claims 
from the theorem are manifestly satisfied following the standard properties of the boson 
operators [Q] 


a' \n) - \n + m), 

( 21 a) 


( 21 b) 


Example 3. The next example will be A in the form of the following disjoint sum: 

A —aQa2a2a4 +a[a^a(,a-y. (22) 

Assume the ground state of A to be = |n)o |n)i 10 ) 2 - 7 . 

AV'-°^ ^ V7i(|n- l)o|n)i |l)f 3^4 |0)f67-I-|n)o |n - l)i [ 0)234 |l)f67 ) = 

We can go on and derive a general expression for and 


Since the state on the right of the previous equation is normalized, the contains a 
square root and A^ could be quite difficult to find (and unnecessary since we need A^/ip 
in the end). So in the next, morally similar, example we derive it using our insight from 
Eq. (fTOl) . 


Example 4. Let 

A-a(ja2 + a{a3 

and = |n)o |n )3 [ 0)2 [ 0 ) 3 . We find 


A'^ ^ jvac) = 

p-i 

m=0 


n 

p — 1 — m 


\n-m)o\n-p + l + m)^\m )2 Ip - 1 - m)^ 

(23) 


By matching the coefficients of the corresponding basis vectors, the calculation 
AA''t/^^^'“^^ ^ = p( 2 n -p + 

reveals the constant of proportionality. We can verify it by a direct calculation for p = 2: 

= (0^02 + a|a 3 )(aoa 2 + aia3)(Vn |n)o jn - l)i [0)2 [1)3 +Vn\n - l)o |n)i [1)2 10)3) 
= ( 4 n- 2 )(v^|n)o|n- 1)1 [0)2 [1)3 +v^|n - l)o |n)i [1)210)3) 

= p(2n-p-Ll)|p^2V'^^^- 

Convinced about the validity of the proposition we can set out to show its main 
application. The sum (A' -l-A)^ contains 2^ summands but when acting on xp, some of 
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them become zero. What is the exact condition for a general summand of +A)^ to 
survive its encounter with the ground state? To find out, we write 


C) 

(A'i' + A)^ = ^ ^ (A )'"!... (Ar^-^(A^r^, (24) 

m=01S|=l 


where 







(25) 


At least one half of all the products on the RHS of Eq. (I24D vanishes upon acting on 
xp because the number of A operators outnumber their conjugates. This is a necessary 
but not sufficient condition to eliminate all the products that do not contribute. For 
example, if there is just one A operator that is, however, the first one to encounter xp, 
the whole summand disappears. A general product in Eq. (I24D is a bosonic word of the 
length k according to (I25D 


(A)'"i(Ai')'"2... (A)'"^-i(A'i')'"2A 


(26) 


From the way it acts on xp it can be seen that the word yields zero if and only if the 
number m 2 i_i of annihilation operators A on the (2i — l)-st position is greater than the 
number m 2 j of all creation operators A' to its right valid for any 1 < i < J. But 
this is precisely the defining property of a Dyck path 2)(k, 0, 62 ). We associate the step 
operators with the lattice steps 

A^D = (1,-1), (27a) 

A'i—[/ = (!, 1) (27b) 


and a bosonic word is nothing else than a Dyck word of instructions of how a Dyck path 
is generated. Indeed, the constraint on a Dyck path to remain confined in the positive 
quadrant is equivalent to the condition for a bosonic word (product) in Eq. 024D to be 
nonzero. The converse is true as well because of the way the annihilation operator acts 
on a ground state and so we have a bijection between the set of Dyck paths Ti{k, 0, 82 ) 
of cardinality G(k, 0, 82 ) and the number of products on the RHS of 

k G(fc,0,52) 

liA^+Afxp= Y, S iArKAT^...IiA)’^^-^liA^r^xp, (28) 

m=[fc/2J 1S|=1 


where 


j j 


II 

M 

1 

M 

3 

to 

1 

(29a) 

j=l 

J=1 


J 

J 



(29b) 

j=l 

J=1 



The first line is extracted from the definition of S (the length of the boson word) and 
the second line is the difference between the number of A and A' determining the y 
component of the end point of the corresponding Dyck path. Note the lower limit of the 
outer sum in Eq. (I28D . It removes the words with a number of annihilation operators 
greater than the number of creation operators. They do not correspond to any Dyck 
path. 
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(b) 


Figure 3. Two examples of the highest (solid) and lowest (dotted) Dyck path for 
a given k, 82 and 5^ = 0. Two different colors distinguish the ascending (black) 
and descending (red) sections of the paths. The letters U; and d; form a Dyck 
word and their role is explained in Sec. 12.31 


We may associate the p-th horizontal line of the integer lattice with the state 
The X axis is then the ground state 'ip. With this insight in hand, we generalize the 
above construction for an arbitrary starting state xp^^^\ where we set 5^ = p: 

k G{k,5i,52) 

(A'^+A)''i/i<^^i^ = ^ ^ (A)'"i(A'^)'”2...(A)'"2J-i(A'i')'"2JT^(:5i)_ (- 30 ) 

m=Lfc/2J-L5iJ |S|=1 

Note that Sj can be greater than k in which case the lower limit of the outer sum is 
zero. 

To illustrate our technique, we will study in this paper the ground state case 5^ = 0 
corresponding to Eq. (I28D . The analysis for an arbitrary ladder state can be easily 
generalized following the approach developed in the next section. 

2.3. Evaluating Dyck paths. In this section we present our main result. Eor that 
purpose we introduce a different kind of Dyck path description equivalent to a Dyck 
word W. Previously we introduced the highest and lowest Dyck paths, D(fc, Sj, ^ 2 ) and 
D(k, 5 ^, ^ 2 ), respectively. They are illustrated in Eig.[^for k = 14, 5 ^ = 0 and 82 = 0, 2 . 
The letter U; denotes the ascending path segment lying between the (i — l)-st and i-th 
horizontal line of the lattice. Similarly, d; denotes the descending segment connecting 
the i-th and (i — l)-st horizontal lattice line. At this point, we treat Uj,dj as letters 
of a different (somewhat less economical) kind of Dyck word we will label w which is, 
nonetheless, equivalent to the usual Dyck word W. A Dyck word w is, like W, read from 
the right and we write down the highest Dyck path W from Eig. (a) as an example: 

W = = did 2 d 2 d 4 d^d(,d-yU-jU(,u^U 4 U 2 U 2 Ui = w. (31) 

Dyck paths are continuous and so two consecutive segments of any Dyck path must be 
one of these four pairs (recall the right-to-left convention): 

Ui+lUi, 
diUi, 

Uidi, 
didi^i- 


(32) 
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It follows that if we have the highest Dyck path, where all Uj are on the right, and there 
is a need to swap the leftmost letter Uj with the rightmost letter dj such that the new 
word represents a Dyck path, the general rule is very simple: 



(33) 


valid for all admissible i. Only then the new word satisfies the conditions in ([32]) no 
matter what letter precedes d; or follows U; on the LHS of (I33D . 

The rule can be immediately put to work by introducing a way of how to generate 
all Dyck paths D(k, 5^, 52 ). We “descend” from the highest Dyck path T){k, 5^, 82 ) to 
the lowest Dyck path ^(k, 0,0) by systematically swapping the highest ascending and 
descending segments U; and d; and generating all Dyck paths ending at (k, 52 ) such 
that no Dyck path is left behind. 

The rule is rather straightforward so let’s illustrate it on the example k = 14 depicted 
in Fig. [ 3 ] (a). The first step is to swap the leftmost U letter of the highest Dyck word 
W = with all the D letters to its left except for the leftmost one. This generates 

the following set of words (the swapped letters are in bold to emphasize the transfor¬ 
mation) : 


D^UDU^ 

D^UDDU^ 

D'^UDDDU^ 

D^UDDDDU^ 

D^UDDDDDU^ 

DUDDDDDDU^. 


(34) 


Looking at the generated Dyck path in Fig. [3] (a) we observe that all Dyck paths reaching 
the level U 5 were obtained. Following the swap rule Eq. (I33D . the same operation using 
the Dyck path representation given by the w word reads 


did2d3d4d5d6U6d6U6U5U4U3U2Ui 

did2d3d4d5U5d5d6U6U5U4U3U2Ui 

did2d3d4U4d4d5d6U6U5U4U3U2Ui 

dld2d3U3d3d4d5dsU(,U5U4U3U2U^ 

did2U2d2d3d4d5d6U6U5U4U3U2Ui 

diUidid2d3d4d5d6U6U5U4U3U2Ui. 


(35) 


Recall that in (I35D we do exactly the same operation as in 034D , just with some more 
information attached to the letters of the Dyck word w. A similar method can be used 
to generate all Dyck paths D(k, 0, 82 ) by starting from T>(k, 0, ^ 2 ). 

In the next step we would have swapped two U’s and let both gradually propagate 
to the last allowed word. But already in this case the resulting words get complicated 
and so we need to simplify our strategy. Before we do so it is perhaps time to reveal 
the reason why. In the crucial step that follows, we associate the letters Uj and dj of all 
Dyck paths, generated according to the rule Eq. (1331) . with the scalar coefficients given 
by Eqs. (|8]) 


di^/J-i 


(36a) 

(36b) 
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resulting in 

i i 

The number x is the result of a transformation we call evaluation of a Dyck path w. 

We already linked boson words Eq. (1261) with Dyck words W via bijection 027D . Be¬ 
cause w is equivalent to W, we automatically linked boson words with Dyck words w. 
But here we go further thanks to the fact that w carries much more information. It 
allows not only to properly count the boson words but also to relate the result of their 
action on a vacuum state ip via the assignment in Eqs. (1361) . 

Now it is clear what our goal is: we collect all w’s, evaluate them and sum them. 
This will be equivalent to the action of a sum of all bosonic monomials with the total 
number of k operators A and A' acting on ip, whose difference between the number of 
creation and annihilation operators is ^2 as a consequence of Eqs. (1361) and (1271) . The 
operations (1361) and 037D are peculiar since we effectively changed the variable’s type 
from a symbol and word to a real number. The number x is not, in general, unique 
since the scalars commute and it may happen that two different Dyck words w 

lead to the same x. But this does not bother us since x will not be used to identify a 
Dyck path. 

In order to execute the transformation Eq. G37D it is critical to be able to list all Dyck 
paths. Here comes the power of Eq. G33D . It informs us that the set of indices i of every 
Dyck path w is identical for the ascending (u;) and descending (dj) section^ Hence to 
systematically list all Dyck paths T){k, 0, 82 ) as the offspring of the highest Dyck path we 
may only focus on the ascending or descending sections of w. But it is very easy to list 
ascending segments only. Every Dyck path 2)(k, 0, 82 ) contains the following ascending 
segments: (see the properties of Dyck paths discussed in Sec. 12.ID . As 

follows from the first row of G32D . the difference of indices of two consecutive ascending 
segments must be one ((i-l-1) —i = 1). This is the greatest possible difference because if 
two ascending segments Uj and Uj are separated by several descending segments {dj.};;. 
it follows from the second and third line of (I32D that j — i < 1. We just rephrased a 
simple fact easily seen by inspecting any Dyck path (see Eigs.[Tl|2]and[^. Now we know 
how to generate all ascending sections of a Dyck path !D(k, 0, 52 ) and we can directly 
substitute A; for U; as dictated by Eq. G36aD and move on to evaluating the whole Dyck 
path. Let’s start with a special case to make the exposition easy to follow. 

Dyck paths D(k,0,0). The product representing the ascending portion of the highest 
Dyck path D(k, 0,0) is The lowest Dyck path corresponds to A^^^. We can 

visualize the products obtained from the ascending segments of all Dyck paths in a 
different kind of diagram in Eig. |4](a) reflecting the example fc = 14 in Eig. [^(a). The 
solid stair path is the highest product and the dotted line is the lowest one. Any product 
from the “wedge” between the two extremal cases comes from an allowed Dyck path iff 
the indices for any consecutive pair A;Aj of a product satisfy 

i-J<l. ( 38 ) 

Using the fact that (i) the number of ascending segments is equal to the number of 
descending segments (and so the number of A products evaluating the ascending sec¬ 
tions is equal to the number of p, products) and (ii) their indices are identical as a 

^Cf. tlie indices of every row in 0351 1 - the indices are the same for u’s and d’s but they are not ordered 
in a particularly lucid way. 
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Figure 4. The x axes count the number of ascending segments. The solid stair¬ 
like path and the dotted line on the left represent the evaluated ascending sec¬ 
tions of the highest Dyck path D(14,0,0) and lowest Dyck path ©(14,0,0), re¬ 
spectively, that are depicted in Fig. [3] (a). In the middle plot we capture the same 
situation for ©(14,0,4) from Fig. [3] (b). In the right plot we reordered the evalu¬ 
ated ascending segments of the lowest Dyck path as suggested by the ascending 
evaluated product for ©(14,0,4) in Eq. 048 D . 


consequence of 033D . the sum over all evaluated Dyck paths reads 


mj+i 


- ni3 

- ^2 

-l - 



^1712-1^^1712-1 



mj=j 


- 7112=2 

- mi=l 

-J - 


We will see in Examples [6l [7] and [ 8 ] and further in Sec. 12.41 that this j-multiple recursive 
sum can be explicitly evaluated (and very fast for large j for that matter). Note that 
the sum’s upper bound does not change: m 2 = = k/2 but for obvious 

reasons we had to introduce distinct dummy indices. Moreover, the number of sums in 
the recursive formula is also j = k/2. So indeed, the outermost sum is not a sum at all 
since visualized by the overlapping segment of the 

solid and dotted line in Eig. |4| (a). Perhaps it is time for an example. 

Example 5. Let k = 8 and this case can be illustrated on a staircase diagram like the 
one in Eig. |4] (a). In the present case the diagram ends with the stair where A 4 lies. 
We can recreate Eq. 039D from “inside” by following the condition g38D for A; Aj and the 
items (i) and (ii) above. Therefore, the innermost sum reads A^/Ti -F A 2 /t 2 +A 3 P 3 -F A 4 /T 4 
(see the horizontal column in Eig. |4| (a) between A^ on the dotted line and A 4 on the 
stair). The sum is multiplied by the components of the preceding column ({Aj/ij}^^^) 
but it depends on j (as g38D dictates) how many summands of the innermost sum are 
allowed to be multiplied - for Ajp^ the previous sum goes from 1 up to j -F 1 : 


+ A 2 P 2 + ■^ 3/^3 + ■ 74 M 4 ) 
+ 72M2(7iMi + A 2/J2 + ^ 3 ^ 3 ) 

+ 7iPi(AiPi -F A2P2)- 


(40) 
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Moving to the next level, only two values and A 2 JU 2 ) are possible. Hence by 

repeating the exact same rules, the third inner sum reads 

+ ■^2^2 + + -^ 4 ^ 4 ) 

+ ■^2M2'^2M2(‘^iMi + ■^2^2 +-^ 3 ^ 3 ) (41) 

+ + ■^ 2 /^ 2 ) 


together with 


Ai/ii A2P2 (^iMi + ^2M2 + ^ 3 /^ 3 ) 

+ AiPiAiPi(AiPi + A2/T2) 

obtained from the second and third line of 04OD by multiplying by A^p^. Finally, the last 
sum (not being a sum since the only possibility is AijUj) terminates the recursion: 

Ai|UiA2JU2^3M3(^iMi + ^2^2 + ^3^3 + ^ 4 ^ 4 ) 

+ AiPi A2/T2A2/i2(‘^lMl + ^ 2^^2 + ■^ 3 /^ 3 ) 

+ AiPiA2jU2^iMi(^iMi + ^ 2 ^ 2 ) (43) 

+ Ai/riAi/riA 2 /T 2 ('^iMi + ■^2^2 + ■^3/^3) 

+ AijUiAijUiAi/ii(AijUi + A2JU2). 


By revising the number of summands (equal to 14), which counts the number of Dyck 
paths D(fc, 0,0), we may verify that it agrees with Eq. © for k = 8 . 


Dyck paths D(fc, 0, 52 ). As already mentioned, the strategy to list all possible Dyck paths 
D(k, 0, 52 ) is similar to the case 52 = 0. Starting from the highest Dyck path we swap 
our way to the lowest Dyck path. Instead of doing it explicitly and encountering con¬ 
fusing expressions like those in Eqs. g34D and (I35D . we better use the rule Eq. (I33D . But 
this time the situation seems to be more complicated. The descending section of any 
Dyck path is shorter than the ascending one whenever 82 > 0. So how do we know 
what Uj and d; belong to the same Dyck path to properly evaluate the path according 
to Eq. (I37D ? Because D(fc, 0, 52 ) is confined to the positive quadrant, the descending 
section of the lowest Dyck path D(k, 0, 82 ) contains (k — 82^/2 of segments 

di. di, (44) 

k-S2 

2 

whose evaluation (following Eq. (I37D ) equals the product 

k-52 

Ml^ . (45) 

The highest Dyck path 2)(k, 0, 52 ) is always of the form 

dg^-i-i. dkj^ . Ui (46) 

2 2 

^-V--V-^ 

^-^2 fe+52 

2 2 

and the number of d; segments correctly coincides with 044D . When (I46D is evaluated, 
we get the expression 

k-52 k+52 

2 2 

n ^52+i n ^47) 

i=l i=l 

As previously described, to generate all Dyck paths we apply the rule Eq. ([33]) on the 
highest Dyck path Eq. (1461 ) and end up with the lowest one (Eq. 044D ). But the rule in 
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Eq. (I33D treats the ascending and descending segments (more precisely, the indices of 
their letters U; and dj) in the same way. So it treats the indices of A; and jU; in the same 
way as well. As a consequence, only the underbraced portion of the ascending segment 
product taken from Eq. 047D 

k+S 2 

2 

I A; = X ... X Ag^ ■^52+1 X • • • X ■^(fc+52)/2 

i=l -•>-- 

k-52 

2 


will change as we descend to the lowest Dyck path. After the lowest Dyck path is 
evaluated, its value must necessarily be equal to the product 

k-52 k-52 ^2 

Ml' ^i' r\^i C48) 

(I) (ii) 

(III) 

The [JLi product denoted (I) is the evaluated descending section, product (II) corre¬ 
sponds to the part of the ascending section whose number (fc — 82)12 must be equal 
to product (I) and product (III) is the rest of the evaluated ascending section. We can 
pictorially see what is happening in Eigs. |4](bi) and (b 2 ) on the example of (D(14,0,4) 
from Eig. [3](b). The solid and dotted staircase in the middle plot is the evaluated as¬ 
cending segment of the highest and lowest Dyck path as can be directly deduced from 
Fig-E](b) using Eq. (I36al) . Similarly to !D(k, 0,0) exemplified in the left panel, all al¬ 
lowed evaluated ascending sections lie between the two extremal cases. An ascending 
section is allowed iff the indices satisfy i — j < 1 for any consecutive pair AjA^ of a 
product. It is, however, extremely advantageous to reorder the lowest evaluated sec¬ 
tion according to the prescription derived in Eq. (I48D . Then we are able to see the 
product of (II) and (III) from Eq. (1481) in the right panel of Eig. |4] as the dotted path. 
More importantly, all products contain expression (III) and so it can be factored out. In 
the example, this happens to be A 1 A 2 A 3 A 4 because 52 = 4. The reason we underwent 
this laborious transformatior@ is to find out the summing formula over all evaluated 
Dyck paths 2)(k, 0, 52 ). With the help of Eq. 048 D and following the strategy leading to 
Eq. G39D we are able to read off the summing formula from Eig. |4](b2): 


! = 1 


''j+l 


I i I 'y I 


■J+1 


■ mj=] 


r m 3 




^2=2 


m 2 


A„ 


-iM; 


mi 


*- 771i=l 


(49) 


Eq. 049 D is a modification of Eq. 039D in a transparent way: the total number of sums is 
j = {k — 82)12 - this is where the solid and dotted paths in Eig. |3] (b 2 ) do not overlap. 
The overall multiplicative constant is the overlapping portion of the diagram where the 
values of A; and jUj are fixed. The upper bounds in 049 D have again the same value 
m 2 = ... = ^ (k-f 82)12. 

Eqs. 039D and 049 D can be used to calculate Eq. 028D . This is where the true power 
of our method is exposed: summing over all evaluated Dyck paths is greatly facilitated 
whenever the action of A and A' in Eqs. dS]) is sufficiently “nice”. But polynomials are 


^The transformation from (bj) to (b 2 ) can also be understood as a consequence of being scalars 
and therefore commuting. Then, in the right panel, the allowed evaluated ascending sections are suitably 
reordered allowed sections from the middle panel. 
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always summable (via Bernoulli’s formula) and so the iterative sums in Eqs. G39D and 
049D are all closed expressions. 

To compare our algorithm with other analytical approaches it is fair to say that for a 
very special case of A = a‘^ the expression J\f [(a^^+a‘^)*^] can be used for the same task. 
The methods to normal order the expression were developed in 11171 IlSIl . A different 
kind of combinatorial insight was used by the authors (more in the spirit of 11161] ) with 
the help of generalized Stirling numbers. The results were subsequently used for the 
calculation of the expectation value {(a^^ + a‘^)^)\a} for a coherent state |a). From the 
conceptual point of view, normally ordered expressions are satisfactory since they are 
essentially closed and hold in the strong sense, that is, without acting on any state. 
But that may be actually a mild disadvantage because it still has to act on a state and 
the resulting formula can become further complicated. What is even perhaps less clear 
is the computational complexity of how effectively it can be evaluated for k » 1 in 
actual calculations. But irrespective of that, here we aim at multiboson operators A 
and A' (constituting physically interesting Hamiltonians), where to the author’s best 
knowledge, analytical results are essentially non-existent. This will be further used to 
evaluate the action of the isometry Eq. ([1]) expanded as a Taylor series to high order in 


k. 

The only (partial) exception seems to be [|2^ where, however, two boson modes are 
coupled such that [ 0 ^, 02 ] = 1. This is not the behavior of ordinary bosons considered 
here and in quantum field theory or quantum optics. 


2.4. Complexity estimation. Let’s analyze the complexity of calculating Eqs. ([39] ) and 
049D for the boson operator A being a product of d creation and annihilation operators. 
This will give us an idea of how tractable our method is. We will focus on 039D as the 
worst case scenario and denote the sums’ upper bounds m. From the general action of 
a boson creation and annihilation operator in Eqs. G21D we find that is a 

polynomial of order d (written symbolically). Then the innermost sum of 039D 

m m 

mi=l mi=l 


is a polynomial of order d -I- 1 and so the cost of calculating the sum is 0(dm). To 
evaluate the following sum, the polynomial is multiplied by and this is 

again a polynomial of order d: — f^'^\ The polynomial is summed over m 2 but 

only to m — 1. In Eq. 039D it is the lower bound that changes by one so it is the same 
number of operations. The whole recursive formula can then be symbolically depicted 
as follows 


m 

S 

fW y(d-Ll) 


m—1 

S 

yr(d-Ll)y(d) ^ j:(2d+2) , 


m times. 


m—2 

S 

jr[2d+2)j:(d] j:(3d+3) 


(51) 
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We first find the order of the resulting polynomial. From the RHS of flSlD we deduce 
that the polynomial after the last iteration m = fc/2 is of order m(d + 1) = k/ 2 {d + 1) 
(recall that since 82 = 0, k is always even). For the complexity estimation we observe 
that the i-th recursive sum on the LHS of flSlD is a sum over m — i + 1 summands and a 
polynomial of order d i + i — 1. Henceforth 

m j 

^(m — i + l)(di + i — 1) = — (^(d + l)m^ + 3dm^ + 2dm — (52) 

i=l 

and so the recursive calculation roughly scales as 0(k^). 

This is the main advantage of our method. It provides a tremendous (exponential) 
speed-up for calculating There are 2^ summands after expanding (A' ±A)^ 

and so their number grows exponentially with k. This is prohibitively slow to compute 
even for a small k and it is perhaps obvious that a similar obstacle remains even after 
we take into account the truncated sum in Eq. (I28D . In more detail, one finds 

Y, G(/c, 0,52)= r fc ] - 0 ( 2 '^-'°®^'^), (53) 

52=0,2,or VI 2 JY 

52=1,3,.. 

where the Stirling formula for the factorial can be used to get the asymptotic estimate 
of the RHS. This is the most favorable case since for 5^ > 0 the number of summands 
in Eq. (1301) is always bigger for fixed k and 82 - By setting 81 = k we get the worst-case 
scenario where the grotvth is simply (9(2^) since none of the summands from (A' -F A)^ 
acting on a state will vanish (see the outer sum of Eq. (1301) ). The largest number of 
Dyck paths in the worst case scenario happens when 5^ = 82 for k even as long as 5^ > 
k/2. In this case Eq. (0 tells us that the number of Dyck paths is G(k, 5^, 5^) = i^^/ 2 )- 
The Dyck paths are then confined to a diamond-like shape illustrated in Eig. [5] (a) for 
5i = ^2 = k = 8. In Eig. [2(b) we see in the corresponding staircase diagram that unlike 
the diagram in Eig. |4| (a), the ascending segments are not approaching. Consequently, 
the number of recursive sums (m = k/2 + 1 = 5 in (I50D in this case) in each step is 
constant. Eollowing the estimate as in Eq. (I51D we arrive at 

m ^ 

(m-Fl)^(di-F i — 1) = — ((d + l)m^ + 2dm^ — m(l — d)). (54) 

We see that the calculation is nearly as effective as for in (I52D . 

Let’s go back to the 5^ = 0 case as that will be our source of examples to come. The 
number of sets of Dyck paths with the same value of 82 grows linearly with k in the 
expression (A' -F A)^. Because we have just shown that the sum over these “52-sets” in 
Eq. (I39D goes as Oik^), the calculation of (1281) becomes tractable. But this also means 
that the calculation of 

K 

= exp [r(A' — A)] ~ ^ 

k=o 

is polynomial in K is the Minkowski vacuum |0)). The number of summands 

in (I55D grows linearly with K but as will be discussed in Example [2 for k = 100, even 
this can be significantly reduced. 

Let’s demonstrate the power of Eqs. 039D and 049D . 
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(a) (b) 

Figure 5. All Dyck paths I>(8,8,8) are confined to the diamond on the left where 
the ascending segments of the highest and lowest (dotted) Dyck path are evalu¬ 
ated by A;. In panel (b) we see the corresponding staircase diagram to illustrate 
the complexity of calculating in the worst case scenario that occurs 

for a given k when 5^ = 82 and 5^ > k/2 {k even). Here we set 5^ = k. 


Example 6. As a generic single-mode example we will calculate (A' +A)^ |0), where 
A = a^, for several values of k. 


k-2 


In this case only two possibilities exist: ^2 = 0,2. For 82 = 0 the number 
summands j in Eq. (I39D equals k/2 = 1 and so there is nothing to sum - only 
the outermost sum survives. But this is not a sum at all as previously discussed 
because there is only one Dyck path D(2,0,0). It remains to calculate (jLi and 
Ai with the help of Eqs. (1211) . Hence, by denoting — |3(k — 1)) we get 

from Eqs. {|8]) 


Afc = jUfc = V3fc(3k-l)(3fc-2) 


(56) 


and so = '/V.. 

For 82 = 2 the situation is similar. The number of sums in Eq. (I49D is j = 
(k — 82)12 = 0 so only the overall factor A 1 A 2 is present. From Eq. (1561) we get 
A 2 = V4 X 5 X 6. The final answer is 


(A^ +Af |0) - Ai/ii |0) +A 1 A 2 |6) = 3! |0) +V^|6). 


k-3 


Q •’ 3 2 

This can be readily confirmed by directly calculating (a^ -Fa' ) |0). 

For 82 = there are two Dyck paths 25(3,0,2). Let’s evaluate this case only. 
The number of summands in Eq. 049 D is now one ((k — 82)12 = 1) and its upper 
bound is m 2 = (k + 82)12 = 2. But to start showing the potential of the Dyck 
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path approach for general k we will not set m 2 = 2 yet. Instead, we calculate 


fc = 100 


m2 


k=l 


l)(3fc - 2) = -m 2 (m 2 + l)(3m2 - 2)(3m2 + 1) 


(57) 


and take the liberty of setting m 2 
final answer given by Eq. 049 D : 


2 in the evaluated sum after we write the 


m 2 (m 2 + l)(3m2 - 2)(3m2 + 1 ) = V^(3! + 4x5x6)-v^ 126. 

This is the first step to achieve a closed expression even for a large k to be used 
in the next two examples. The calculation again agrees with the |3) coefficient 
of 

(a^'^ + |0) = y^(3! + 120) |3) +\f9\ |9). 

Here we will illustrate the appearance of one nested sum in Eq. 049D for ^2 = 1 
coming from j = 2. The recursive summation results in 

9 


1120 
X (2835mf^ 


(m3 - I)m3(m3 + l)(m3 + 2) 


5130mo 


3915m3 + 8034m3 


808), 


(58) 


L3 

where now we set m 3 = 3. The final answer for 62 = 1 is therefore 76356. 
This can be verified by hand but it is already rather lengthy. 

To show the strength of our method we calculate |0) = (a'^ + 

a^) |0) for 82 = 0. The number of Dyck paths given by Eq. © is 

1978 261657 756160 653 623 774456. 


fc = 200 


The brute-force calculation is out of the question on today’s computers. 

Using our approach we see that Eq. (I39D contains 50 recursive sums. It takes 
roughly 20 seconds to obtain an analytical result in Mathematica on a single¬ 
core average laptop. The result (the equivalent of (1581) ) is a polynomial of 
order 200 whose printout would take 24 A4 pages. The order coincides with 
the derivation below Eq. (1511) for d = 3 (note k/2(d + 1) = 200). Given the re¬ 
cursive character of Eq. (I39D . the code used to generate the result is a one-liner. 
To calculate the complete expression (a' -I- a^) |0) we, of course, have to 

calculate the remaining cases for the rest of the 52 ’s (62 = 2 ,4,..., 100 ) given 
by Eq. 049 D . The number of recursions is j = (k — 62)72 but these recursions 
are already contained in the largest calculation for 62 = 0. The only thing 
that differs for each 62 is the sums’ upper bound. Recall that we calculate the 
sum for an arbitrary upper bound m^, where 2 < i < j + 1 , and plug the value 
(k -I- 62)72 after the calculation. So we only need to save all j intermediate 
nested sums calculated for 62 = 0 , plug mj = (k + 62)72 for the rest of 62 and 

multiply it by the overall factor Ofli required by Eq. 049 D . Indeed, the 
first nested sum in the k = 100 calculation is precisely Eq. 058D . where instead 
of m 3 we have m^Q. So all the remaining 62 operations are computationally for 
free! 

This case is to verily the cubic complexity scaling (for a fixed d) derived in 
Eq. 052D for the recursive sum Eq. 039D . It takes some 240 seconds to obtain an 
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analytical result in Mathematica. This is indeed expected: 200^/100^ = 8. The 
result is a polynomial of order 400 whose printout fits 105 A4 pages. 

Even without using the computational shortcut for ^2 > 0 described in the fc = 100 
case, we already got rid of the exponential time complexity by summing over all Dyck 
paths T>{k, 0 , 82 )- 

The cubic boson expression from the previous examples can be thought as a Hamil¬ 
tonian generating a unitary evolution operator, see Example 0for the exact formulation. 
But the following unitary is even more important: the two-mode squeezing operator. 
Its factorization is a routine procedure and so we can compare it with our approach. 

Example 7. Let’s investigate — exp[r(a|a 2 “ where is a ground 

state coinciding with the Minkowski vacuum |0) annihilated by U;. We will compare the 
closed form of the amplitude coefficients provided by 

2 00 

exp[r(a|a 2 - 0102 )] |0) ^ —— Vtanh” r |n)i |n )2 (59) 

cosh r 

n=0 

with 

K j.k 

^|0) =^^—(aia2“ “ 1 ^ 2 / 10 ) (60) 

fc=o 

coming from Eq. (1551) . The calculation reduces to (aju^ ~ |0) = (A —A')^ |0) for 

0 <k< K. Using the methods developed in this paper we can afford to set K very high 
and study the evolution for any fixed value of r for which V converges. Erom Eq. (1211) 
we find 

M/c ^ 

by setting = \k — l)i \k — 1)2 in Eqs. d^. It is perhaps clear that for any k, 

the admissible values of 82 tell us to what state \n)i\n )2 we calculate the amplitude 
contribution. So, for example, whenever 82 = 0, the contribution goes to the vacuum 
amplitude |0), the ^2 = 1 contributions go to |1)^ | 1)2 etc., up to \k)i \k )2 since we know 
that 0 < ^2 < k if k is even (if it is odd we offset ^2 by one, see Sec. 12.ID . 

If we set K = 200, all 200 calculations (aju^ ~ ^ 1 ^ 2 )^ |0) t^ke roughly 250 seconds to 
complete without using further optimization described in the k = 100 case of Example!^ 
As a matter of fact, the result of the k-th expansion summand could have been used to 
speed up the calculation in all the exponents that followed. Neither this optimization 
was implemented. Just for comparison, considering only the highest contribution (k = 
200) to the vacuum amplitude (^2 = 0), they come from the immense number 

r m 

G200 ~ ^ 

of Dyck paths, according to Eq. (0. 

If the evolution operator V was evaluated using its Taylor expansion, it would be 
necessary to assume r —> 0 as only a few first expansion coefficients are possible to 
manage. Not in our case. Eor K = 200 and a very higf0 value of the coupling constant 

^The magnitude of r is a relative notion especially because 0 < r < 00 . Nevertheless, compared to what 
has been possible so far for the Taylor expansion, it is indeed high and following our performance analysis 
in Sec. 12.41 we can now calculate any finite value just by sufficiently increasing the Taylor expansion. From 
a physical point of view it is also high. If translated to the quantum optical scenario, where r is a squeezing 
parameter, then the squeezing value r = 1 corresponds to the noise reduction lOlogjge^'' 8.7 dB. This 
is roughly current state-of-the-art in optical experiments ifi^] . If we use a special-relativistic comparison 
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r — 1, we get very good agreement with tanh'’ r / cosh r up to n = 33: 
tanh^^ 1 

033 --;-~ 0.000081001 vs. 033 = 0.0000799909. (61) 

coshl 

The lower coefficients starting with n = 30 are essentially indistinguishable from the 
ones given by the analytical expression. By comparing the Taylor expansion of the an¬ 
alytical coefficients o^ = we indeed see the perfect match. Also recall 

that the calculation of Eq. (1601) is perturbative but analytical and the value of the cou¬ 
pling constant r is inserted after the calculation. So the expansion is not needed to be 
recalculated for every numerical value of r. 


Example 8. The evolution operator 

W = exp [r(a'^ — a^)] 


(62) 


is supposed to describe the three-photon degenerate parametric amplifier. Some time 
ago it attracted a lot of attention by an observation made in [143n that the vacuum 
expectation value does not converge for r > 0. The situation was clarified by Braunstein 
and McLachlan who uncovered the immediate reason for divergence and used an 
analytic continuation method based on Fade approximants showing convergence 
for a finite interval of r. It is possible that the cause of the bad behavior of W is the 
lack of a self-adjoint extension on the domain of interest of the Hamiltonian = 
r(a' — a J as this is a subtle issue for unbounded operators. But to the author’s 
knowledge this problem has not been clarified yet. What is clear is that the Taylor series 
is essentially useless but even here our method can be used. Eollowing the footsteps 
of llsbll we converted the Taylor expansion coefficients of is the Minkowski 

vacuum |0) and (1561) was used in the calculation and the difference in the minus sign 
was taken into account as well) to the Fade polynomial [83/ 831 and arrived at an 
identical result. See Fig. and compare it with Fig. 2(a) in [136I1 . 


Example 9. For our final example of the developed method I invite the reader to 
study [l45n , where a recently developed model of a unitary black hole evaporation based 
on the trilinear boson Hamiltonian 


Htri = r(ab'c'— a'be) (63) 

is investigated. Similarly to the expression in Eq. (1621) the unitary operator W = 
exp [Htril is not known to be easily factorizable but unlike the cubic evolution oper¬ 
ator it exhibits no problems when it comes to the convergence of its Taylor series. 


3. Conclusions 

In this work, we introduced a method of evaluating evolution operators for the inter¬ 
action picture Hamiltonians of the form H = r(A' — A), where A is a multimode boson 
monomial and an extensive class of sums of boson monomials and r is a coupling con¬ 
stant. The calculation of the evolution operator V = exp [r(A' — A)] is in general a 
difficult problem. Even though a factorization is always possible, for example, by virtue 
of the Zassenhaus formula or other decoupling techniques, a simple factorization into 
a small number of products is available only if the operators A and A' satisfy favorable 

based on the local isomorphism St/(1,1)7^2 SO(2,1) lll^ldil] . then 2r = 2 is the rapidity corresponding 

approximately to 96.4% of the speed of light. 
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Figure 6. The vacuum expectation value for the cubic boson evolution operator 
W = exp [r(a^^ — a^)] is plotted as a function of the coupling constant. It coin¬ 
cides with the result of [13611 showing that our method can be useful even for the 
evolution operators whose Taylor series diverges. 


algebraic properties, such as vanishing commutators. If this is not satisfied, it may hap¬ 
pen that there is no simple method left and the remaining techniques at our disposal 
are computationally demanding. 

The Taylor expansion in r with its subsequent action on a state of interest (most 
often a ground state of A defined = 0) is certainly not the most efficient method 

as the number of summands increases exponentially with the expansion order. The 
method developed here overcomes this problem and allows us to analytically calculate 
the action of the evolution operator V to a very high order of its Taylor expansion and 
previously hardly accessible values of the coupling constant. This is possible due to 
an insight that a combinatorial structure known as a Dyck path can be related to the 
action of boson monomials on a ground state tp. It helps us to dramatically reduce the 
complexity of analytically calculating (A' from O to 0{dk^) for all 

multimode boson monomials A of the length d and any ground state The boson 
expression is therefore not required to satisfy any special algebraic property except 
for the fact that they behave as ladder operators: oc 'ip^P\ The state xp^P^ is 

any state from the tower of states generated by the repeated action of A' on a ground 
state. The result is applicable for V acting on a ground state but the technique is 
equally efficient for all 'ip^P^. The speed-up is possible due to the existence of a summing 
recursive formula that can be explicitly evaluated for all such A and A' forming (A' ±A)^. 
Consequently, the complexity of calculating the evolution operator V'lp^^^ is polynomial¬ 
time as well. All boson evolution operators (their Taylor expansions) are expressed in 
terms of A and A' and we believe that our technique could be relevant for computational 
purposes in condensed matter, quantum field theory, quantum optics and other branches 
of quantum physics. We illustrated the method on a non-trivial example of A and also 
verified it in the case of a two-mode squeezed operator, where a simple factorization 
procedure is well knovra and the degenerate cubic Hamiltonian, whose Taylor series is 
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known to diverge for any value of the coupling constant r and has to be analytically 

continued using, for instance, the Fade approximants. 
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